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Abstract. The steady state of a dilute gas enclosed between two infinite par- 
allel plates in relative motion and under the action of a uniform body force 
parallel to the plates is considered. The Bhatnagar-Gross-Krook model ki- 
netic equation is analytically solved for this Couette— Poiseuille flow to first 
order in the force and for arbitrary values of the Knudsen number associated 
with the shear rate. This allows us to investigate the influence of the exter- 
nal force on the non-Newtonian properties of the Couette flow. Moreover, the 
Couette— Poiseuille flow is analyzed when the shear-rate Knudsen number and 
the scaled force are of the same order and terms up to second order are retained. 
In this way, the transition from the bimodal temperature profile characteristic 
of the pure force-driven Poiseuille flow to the parabolic profile characteristic 
of the pure Couette flow through several intermediate stages in the Couette— 
Poiseuille flow are described. A critical comparison with the Navier— Stokes 
solution of the problem is carried out. 

1. Introduction. Two paradigmatic stationary noncquilibrium flows are the plane 
Couette flow and the Poiseuille flow. In the plane Couette flow the fluid (henceforth 
assumed to be a dilute gas) is enclosed between two infinite parallel plates in relative 
motion, as sketched in Fig. 1(a). The walls can be kept at different or equal tem- 
peratures but, even if both wall temperatures are the same, viscous heating induces 
a temperature gradient in the steady state. If the Knudsen number associated with 
the shear rate is small enough the Navier-Stokes (NS) equations provide a satis- 
factory description of the Couette flow. On the other hand, as shearing increases, 
non-Newtonian effects (shear thinning and viscometric properties) and deviations 
of Fourier's law (generalized thermal conductivity and streamwise heat flux com- 
ponent) become clearly apparent [20]. These nonlinear effects have been derived 
from the Boltzmann equation for Maxwell molecules [14, 28, 34, 40, 57], from the 
Bhatnagar-Gross-Krook (BGK) kinetic model [6, 19, 43], and also from general- 
ized hydrodynamic theories [49, 51]. A good agreement with computer simulations 
[16, 17, 25, 30, 31, 37] has been found. The plane Couette flow has also been ana- 
lyzed in the context of granular gases [60, 65]. In the case of plates at rest but kept 
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Figure 1. Sketch of (a) the Couette flow, (b) the force-driven 
Poiscuillc flow, and (c) the Couette-Poiseuille flow. 



at different temperatures, the Couette flow becomes the familiar plane Fourier flow, 
which also presents interesting properties by itself [4, 16, 17, 24, 29, 33, 40, 41, 42]. 

The Poiseuille flow, where a gas is enclosed in a channel or slab and fluid motion is 
induced by a longitudinal pressure gradient, is a classical problem in kinetic theory 
[10, 35]. Essentially the same type of flow field is generated when the pressure 
gradient is replaced by the action of a uniform longitudinal body force F = mgx 
(e.g., gravity), as illustrated in Fig. 1(b). This force-driven Poiseuille flow has 
received a lot of attention both from theoretical [2, 3, 15, 18, 21, 32, 38, 39, 44, 
49, 51, 55, 56, 58, 59, 64, 67] and computational [1, 12, 22, 23, 27, 38, 61, 62, 
68] points of view. This interest has been mainly motivated by the fact that the 
force-driven Poiseuille flow provides a nice example illustrating the limitations of 
the NS description in the bulk domain (i.e., far away from the boundary layers). 
In particular, while the NS equations predict a temperature profile with a flat 
maximum at the center, computer simulations [27] and kinetic theory calculations 
[55, 56] show that it actually has a local minimum at that point. 

Obviously, the Couette and Poiseuille flows can be combined to become the 
Couette-Poiseuille (or Poiseuille-Couette) flow [9, 36, 50, 53]. To the best of our 
knowledge, all the studies on the Couette-Poiseuille flow assume that the Poiseuille 
part is driven by a pressure gradient, not by an external force. This paper intends 
to fill this gap by considering the steady state of a dilute gas enclosed between two 
infinite parallel plates in relative motion, the particles of the gas being subject to 
the action of a uniform body force. This Couette-Poiseuille flow is sketched in Fig. 
1(c). We will study the problem by the tools of kinetic theory by solving the BGK 
model for Maxwell molecules. 

The aim of this work is two-fold. First, we want to investigate how the fully 
developed non-Newtonian Couette flow is distorted by the action of the external 
force. To that end we will assume a finite value of the Knudscn number related to 
the shear rate and perform a perturbation expansion to first order in the force. As a 
second objective, we will study how the non-Newtonian force-driven Poiseuille flow 
is modified by the shearing. This is done by assuming that the shear-rate Knudsen 
number and the scaled force are of the same order and neglecting terms of third 
and higher order. In both cases we are interested in the physical properties in the 
central bulk region of the slab, outside the influence of the boundary layers. 
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The organization of the paper is as follows. The Boltzmann equation for the 
Couette-Poiseuille flow is presented in Sec. 2. Section 3 deals with the NS descrip- 
tion of the problem. The main part of the paper is contained in Sec. 4, where the 
kinetic theory approach is worked out. Some technical calculations are relegated to 
Appendix A. The results are graphically presented and discussed in Sec. 5. The 
paper ends with some concluding remarks in Sec. 6. 

2. The Couette Poiseuille flow. Symmetry properties. Let us consider a 
dilute monatomic gas enclosed between two infinite parallel plates located at y = 
±L/2. The plates are in relative motion with velocities U± along the x axis and 
are kept at a common temperature T w . The imposed shear rate is therefore u = 
(U+ — U-)/L. Besides, an external body force F = mgx, where m is the mass of a 
particle and g is a constant acceleration, is applied. The geometry of the problem 
is sketched in Fig. 1(c). In the absence of the external force (g = 0) this problem 
reduces to the plane Couette flow [see Fig. 1(a)]. On the other hand, if the plates 
are at rest (ui = 0), one is dealing with the force-driven Poiseuille flow [see Fig. 
1(b)]. The general problem with w/0 and g ^= defines the Couette-Poiseuille 
flow analyzed in this paper. 

In the steady state only gradients along the y axis are present and thus the 
Boltzmann equation becomes 



where / is the one-particle velocity distribution function and J[f, f] is the Boltz- 
mann collision operator [7, 8, 11, 13, 26, 48], whose explicit expression will not be 
written down here. The notation /(?/, v|w, g) emphasizes the fact that, apart from 
its spatial and velocity dependencies, the distribution function depends on the inde- 
pendent external parameters uj and g. As said above, g = and u = correspond 
to the Couette and Poiseuille flows, respectively. 

In general, Eq. (1) must be solved subjected to specific boundary conditions, 
which can be expressed in terms of the kernels K±(v, v') defined as follows. When 
a particle with velocity v' hits the wall at y = L/2, the probability of being reemitted 
with a velocity v within the range dv is K+(v, v')dv; the kernel K-(v, v') represents 
the same but at y = —L/2. The boundary conditions are then [13] 



e(±v y )\v y \f(y = TL/2,v) = Q(±v y ) / dv' K|/Mv, v')e(T<)/(» = TL/2,v'), 



where Q(x) is Hcaviside's step function. In the case of boundary conditions of 
complete accommodation with the walls, so that K±(v,v') = K±(v) does not 
depend on the incoming velocity v', the kernels can be written as 



where tp^ (v) represents the probability distribution of a fictitious gas in contact with 
the system at y = I fL/2. Equation (3) can then be interpreted as meaning that 
when a particle hits a wall, it is instantaneously absorbed and replaced by a particle 
leaving the fictitious bath. Of course, any choice of V 3 t( v ) mus t be consistent with 
the imposed wall velocities and temperatures, i.e. 




(1) 



(2) 



K T (v) = A 




(3) 




(4) 
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k B T w = im J ' dv (v - U T )V T ( V )- (5) 
Inserting Eq. (3) into Eq. (2), one has 

e(±«,)/(y = Ti/2, v) = G(±v y )n T ^ T (v), (6) 

where 

/dVe(T^)|^|/(y = TV2,vQ 

nT - /dv'e(±^)Ki^ T (v') • u 

The boundary conditions (6) are usually referred to as diffuse (or stochastic) bound- 
ary conditions. The simplest and most common choice is that of a Maxwell- 
Boltzmann (MB) distribution [52]: 

3/2 

(8) 



m(v - U T ) 



2k B T u 



The first few moments of / define the densities of conserved quantities (mass, 
momentum, and temperature) and the associated fluxes. More explicitly, 



n(y\uj, g) = J dv f(y, v\u, g), (9) 
n{y\u,g)-u{y\u,g)= \ dvvf(y,v\u,g), (10) 



in 



■n(y\u,g)k B T{y\g) = p(y\u,g) = y / dvV 2 (y, v\u,g)f{y, v\w,g), (11) 
p ij{y\u,g) = m / dvVi(y,v\u,g)Vj(y,v\u),g)f(y,v\u),g), (12) 



q(2/Kff) = -2 / dvV (y,v\u,g)Y(y,v\w,g)f{y,v\w,g). (13) 

In these equations n is the number density, u is the flow velocity, 

Y(y,v\a,u,) = \-u(y\w,g) (14) 

is the peculiar velocity, T is the temperature, k B is the Boltzmann constant, p is 
the hydrostatic pressure, is the pressure tensor, and q is the heat flux. Taking 
velocity moments in both sides of Eq. (1) one gets the following exact balance 
equations 

dyPyy = 0, (15) 

dyPxy = rung, (16) 
d y q y + P xy d y u x = 0. (17) 
Henceforth, without loss of generality, we will assume u x (0) = 0. In other words, 
we will adopt a reference frame solidary with the flow at the midpoint y = 0. 

The symmetry properties of the Couette-Poiscuillc flow imply the following in- 
variancc properties of the velocity distribution function: 

f(y,v x ,v v ,v z \u],g) = f(-y,-v x ,-v y ,v z \uj,-g) 
= f(-y,v x ,-v y ,v z \-w,g) 
= f(y>v x ,v y ,-v z \w,g), (18) 
As a consequence, if xCz/l^iS) denotes a hydrodynamic variable or a flux, one has 

x(y\v,g) = S g x{-y\u,-g) 

= SwX{-v\-u,g), (19) 
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*vv 


+ 






q x 


+ 


% 





Table 1. Parity factors S g and 5 W for the hydro-dynamic fields 
and the fluxes [see Eq. (19)]. 



where S g = ±1 and 5 W = ±1. The parity factors S g and S u for the non-zero 
hydrodynamic fields and fluxes are displayed in Table 1. In general, if % is a moment 
of the form 

X(.v\v,g) = JdvV x k -(y\uj,g)v^v^f(y,v\LJ,g) (20) 

then S g = (-l) k *+ k v and S u = (-l) fc ». 

The general solution to the stationary Boltzmann equation (1) with the boundary 
conditions (6) can be split into two parts [45, 46, 47]: 

/(y,v|w,.g) = f H (y,v\uj,g) + f B {y,v\uj,g). (21) 

Here, fa represents the hydrodynamic, Hilbert-class, or normal contribution to 
the distribution function. This means that fn depends on space only through a 
functional dependence on the hydrodynamic fields, i.e. 

f H (y) = f H [n,T,u x ]. (22) 

The contribution /b represents the boundary-layer correction to fn, so that / = 
fn + /b verifies the specified boundary conditions. The correction fg is appreciably 
different from zero only in a thin layer (the so-called boundary layer or Knudscn 
layer), adjacent to the plates, of thickness of the order of the mean free path. Conse- 
quently, if the separation L between the plates is much larger than the characteristic 
mean free path, there exists a well defined bulk region where the boundary correc- 
tion vanishes and the distribution function is fully given by its hydrodynamic part. 
In the boundary layers the hydrodynamic profiles are much less smooth than in 
the bulk domain. The values of the flow velocity near the walls are different from 
the velocity of the plates (velocity slip phenomenon), i.e. u x (y = ±L/2) / U±. 
Besides, the extrapolation of the velocity profile in the bulk to the boundaries, 
Ux,n{y = ±£/2), is also different from both the actual values u x (y = ±L/2) and 
the wall velocities U±. Of course, an analogous temperature jump effect takes place 
with the temperature profile. The boundary contribution fg for small Knudscn 
numbers has been analyzed elsewhere [8, 47]. 

In the remaining of this paper we will focus on the hydrodynamic part fn (and 
will drop the subscript H), with special emphasis on the corresponding hydrody- 
namic contributions to the momentum and heat fluxes. In order to nondimension- 
alize the problem, we choose quantities evaluated at the central plane y = as 
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units: 

f *l *, *n \ * v lk B T(0) 

f (s,v \a,g ) = -t-—f(y,vuj,g), v = — — , v T (0 = \/ , 23 

» s«>« = 77y— , T ) = ^Tvv; , P{sa,9) = 777—' ( 24 ) 

n(0) T(0) p(0) 

W»,g) = Ko) , q («h* ) = -^^y , (25) 
_ _J_ <9w* * = 9 * = ,„ fi . 

In the above equations we have found it convenient to introduce the dimcnsionlcss 
scaled spatial variable 

S ^ s ^m /V^y), (27) 

«t(0) Jo 

where v{y) is an effective collision frequency. For the sake of concreteness, we choose 
it as 

= (28) 

where rj is the NS shear viscosity. The change from the boundary-imposed shear 
rate ui to the reduced local shear rate a is motivated by our goal of focusing on 
the central bulk region of the system, outside the boundary layers. Note that a 
represents the Knudsen number associated with the velocity gradient at y = 0. 
Likewise, g* measures the strength of the external field on a particle moving with 
the thermal velocity along a distance on the order of the mean free path. 
The relationship (27) can be inverted to yield 

y*^ = f -^frv y* = J/ f nv (29) 

Jo v (« ) v T {0)/v{0) 
The invariance properties (18) translate into 

f*(s,v*,v*,v* z \a,g*) = f*(-s,-v*,-v*,v*\a,-g*) 

{s,v x , - v y, v z \ ~ a,g ) 
= r(-s,v* x ,v* y ,-v* z \a,g*). (30) 

Given the symmetry properties (30), we can restrict ourselves to a > and g* > 
without loss of generality. 

3. Navier Stokes description. To gain some insight into the type of fields one 
can expect in the Couettc-Poiscuillc flow, it is instructive to analyze the solution 
provided by the NS level of description. In the geometry of the problem, the NS 
constitutive equations are 

P XX = Pyy = Pzz = P, (31) 
Pxy = -r)d y u x , (32) 
q x = 0, (33) 
q y = -Kd y T, (34) 

where 77 is the shear viscosity, as said above, and k is the thermal conductivity. 
Inserting the NS approximate relations (31)-(34) into the exact conservation equa- 
tions (15)— (17) one gets 

p = const, (35) 
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(ridy) 2 u x = -rjmng, (36) 

^(r,d y ) 2 T=-(r,d y u x f, (37) 

where Pr = (5/cs/2m)?7/ft ~ | is the Prandtl number. In dimensionless form, Eqs. 
(36) and (37) can be rewritten as 

%ul(8) = -£&g*, (38) 

d 2 s T*(s) = -^ld s u* x ( S )} 2 . (39) 

For simplicity, let us assume that the particles are Maxwell molecules [7, 11, 63], 
so v{y) oc n(y) and v*(s) = n*(s). In that case, Eqs. (38) and (39) allow for an 
explicit solution: 

u x ( S \a,g*) = as-±g*s 2 , (40) 

T*(s\a,g*) = 1 - ^s 2 (6a 2 - 4ag*s + g* 2 s 2 ) . (41) 

Here we have applied the Galilean choice u x (0) = and the symmetry property 
dyT\ y=0 = 0. 

Equation (40) shows that, according to the NS approximation, the velocity field in 
the Couette-Poiscuillc flow is simply the superposition of the (quasi) linear Couette 
profile and the (quasi) parabolic Poiseuille profile. In the case of the temperature 
field, however, apart from the (quasi) parabolic Couette profile and the (quasi) 
quartic Poiseuille profile, a (quasi) cubic coupling term is present. Here we use the 
term "quasi" because the simple polynomial forms in Eqs. (40) and (41) refer to 
the scaled variable s. To go back to the real spatial coordinate y one needs to make 
use of the relationship (27), taking into account that for Maxwell molecules v oc n. 
Instead of expressing s as a function of y it is more convenient to proceed in the 
opposite sense by using Eq. (29). Since 1/v* =T* one simply has 



y*{s) = s 



i f r 2/o2 * , 1 *2 2 

1 s 2a — an s H — q s 

30 V J 5 y 



(42) 



For further use, note that, according to Eq. (41), 



d 2 T* 



dy* 



= -Pr^a 2 . (43) 



Thus, the NS temperature profile presents a maximum at the midpoint y* = 0. 

Before closing this section, let us write the pressure tensor and the heat flux 
profiles provided by the NS description: 

Kx(s\a,g*) = P* y {s\a,g) = PZ,{s\a,g) = 1, (44) 

p xy{ s \a,g*) = -a + g*s, (45) 

q*(s\a,g*)=0, (46) 

q* y (s\a,g*) = s(a 2 -ag*s + ^g* 2 S 2 ^j . (47) 
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4. Kinetic theory description. Perturbation solution. Now we want to get 
the hydrodynamic and flux profiles in the bulk domain of the system from a purely 
kinetic approach, i.e., without assuming a priori the applicability of the NS con- 
stitutive equations. To that end, instead of considering the detailed Boltzmann 
operator J[f,f] we will make use of the celebrated BGK kinetic model [5, 7, 66]. 
In the BGK model Eq. (1) is replaced by 

( d d \ 

V y dy + 9 dv~) -^ y ' v l w ' 5 ) = -K2/Kfl0 [/(J/> V K 9) - M(y, v\uj,g)] , (48) 
where v is the effective collision frequency defined by Eq. (28) and 

\ 3 / 2 / ta2 

m \ ' / mv 



is the local equilibrium distribution function. In terms of the dimcnsionless variables 
introduced in Eqs. (23)-(27), Eq. (48) can be rewritten as 

(l + v* y d s ) r(s,v*\a,g*) = M*(s,v*\a,g*) - f *|o,g*). (5Q) 

x " ' i/*(s\a,g*) ov* 

Its formal solution is 



/*(v*) = (i + v* y d s y 

oo 



fc=0 



a* d 
v* ov* 
q* d 
v ' v* dv* K 1 



(51) 



The formal character of the solution (51) is due to the fact that /* appears on the 
right-hand side explicitly and also implicitly through A4* and v* . The solvability 
(or consistency) conditions are 



dv* 



{l,v>* 2 }r( S ,v*|a, 5 *) = |dv*{l,v>* 2 }>T( S ,v*| a , ff *). (52) 

Let us assume now that g* is a small parameter so the solution to Eq. (50) can 
be expanded as 

f*(s,v*\a,g*) = fZ(s,v*\a) + ft(s, v*\a)g* + / 2 *(s, v» 3 * 2 + • ■ • . (53) 
Likewise, 

X*(s\a, g*) = X *o(s\a) + X*Ma)g* + X^Mff* 2 + ■ ■ • , (54) 
where x* denotes a generic velocity moment of /*. The expansions of n* , u*, and 
T* induce the corresponding expansion of M* . The expansion in powers of g* 
allows the iterative solution of Eq. (51) by a scheme similar to that followed in Ref. 
[54] in the case of an external force normal to the plates. 

4.1. Zeroth order in g. Pure Couette flow. 

4.1.1. Finite shear rates. To zeroth order in g*, Eqs. (50) and (51) become 

(1 + v* y d s ) / *(s,v*|a) = M* (s,v*\a), (55) 

OO 

/o*(s>v*|a) = J2(-<d s ) k M* (s,v*\a), (56) 

fe=0 

where 



(2-kY< a T* 
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These are just the equations corresponding to the pure Couette flow. The complete 
solution has been obtained elsewhere [6, 20, 25] and so here we only quote the final 
results. The hydrodynamic profiles are 

pgHa) = l, (58) 
u* j0 (s|a) = as, (59) 
r *( S | a ) = l- 7 (a) S 2 , (60) 
where the dimcnsionlcss parameter 7(a) is a nonlinear function of the reduced shear 
rate a given implicitly through the equation [6, 20] 

,*a(7)" 
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3 + 2 



where the mathematical functions F r (x) are defined by 
2 



F (x) = \J™ dtte-*/ 2 K (2x- l 'W a ) t F r {x) = (±xj F (x), 



(61) 



(62) 



Kq(x) being the zeroth-ordcr modified Bcsscl function. Equation (62) clearly shows 
that F r (x) has an essential singularity at x — and thus its expansion in powers of 
x, 



F r (x) 



OO 

E 

k=0 



(fc + l) r (2fc + l)!(2/c + l)!!(-a;) A 



(63) 



is asymptotic and not convergent. However, the series representation (63) is Borel 
summable [6, 25], the corresponding integral representation being given by Eq. (62). 
The functions F r (x) with r > 3 can be easily expressed in terms of Fq(x), Fi(x), 
and *2(a;) as 

1-F (x) „ _ 1 



F r (x) 



1 r_3 

— Y 

8x ^ 



F 3 (x) 

r-3 
m 



8x 

(-l) m+r F m (x) 



* 2 (x) - -F r {x), 



(64) 



1 



F r -i(x) — -F r _2(x), r>4. (65) 



It is interesting to compare the hydrodynamic profiles with the results obtained 
from the Boltzmann equation at NS order (see Sec. 3). We observe that Eq. (58) 
agrees with Eq. (35) and Eq. (59) agrees with Eq. (40) for g* = 0. On the other 
hand, Eq. (41) with g* = differs from Eq. (60), except in the limit of small shear 
rates, in which case 7(a) ^a 2 (Note that Pr = 1 in the BGK model). 

The relevant transport coefficients of the steady Couette flow are obtained from 
the pressure tensor and the heat flux. They are highly nonlinear functions of the 
reduced shear rate a given by [6, 19, 20, 30] 

1 + 47^(7) + F 2 (7)], 

P; y , (s\a) = l-2 7 [F 1 ( 7 )+2F 2 ( 7 )], 

P* 0(40 = 1-27*1(7), 



xy,0 



(s\a) = -0*0(7), 



(66) 
(67) 
(68) 
(69) 



ll,o(s\a) 



*o(7) - 1 - 107*1(7) - 87*2(7) 
Qy,o( s \ a ) = a 2 *o(7)s- 



*Mt) 

*i(7) 



(70) 
(71) 



Notice that, although the temperature gradient is only directed along the y axis 
(so that there is a response in this direction through g*), the shear flow induces a 
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nonzero x component of the heat flux [19, 20, 30, 37]. Furthermore, normal stress 
differences (absent at NS order) are present. Equations (69) and (71) can be used to 
identify generalized nonlinear shear viscosity and thermal conductivity coefficients. 

In general, the velocity moments of degree k of /q are polynomial functions of the 
spatial variable s of degree k — 2. An explicit expression for the velocity distribution 
function /q has also been derived [20, 25]. 

4.1.2. Limit of small shear rates. The coefficient 7(a) characterizing the profile of 
the zeroth-ordcr temperature T * is a complicated nonlinear function of the reduced 
shear rate a, as clearly apparent from Eq. (61). Obviously, the zeroth-order pressure 
tensor and heat flux given by Eqs. (66)-(71) inherit this nonlinear character. 

It is illustrative to assume that the reduced shear rate a is small so one can 
express the quantities of interest as the first few terms in a (Chapman-Enskog) 
series expansion. From Eqs. (61)-(71) one obtains 

^4( i+ ! a2+ -)< < ra > 

^o(*H = l + ^(l-^« 2 + -), (73) 

^,oH«) = i-y(i-^ a + ■•■)' ( 74 ) 

P: z , (s\a) = l-^(l-^a 2 + .-.y (75) 

P: yfi (s\a) = -a(l-™a 2 + ---y (76) 

. / , n 14a 3 / 1836 o \ 
<o(*|o) = — ( 1_ rnf + '" ) s ' ( 77 ) 

<4o(s|«H « 2 (l - f « 2 + ■•■)*• (78) 

The terms of order a 2 , a, and a 2 in Eqs. (72), (76), and (78), respectively, agree 
with the corresponding NS expressions, Eqs. (41), (45), and (47). On the other 
hand, as noted above, the normal stress differences (P* x — P* y and P* z — P* y ) and 
the streamwise heat flux component q* reveal non-Newtonian effects of orders a 2 
and a 3 , respectively. 

4.2. First order in g. Couette Poiseuille flow. 

4.2.1. Finite shear rates. To first order in g* Eq. (51) yields 

A*(v*) -M\{v*) = A«(v*) + A( n >(v*), (79) 

where 

OO OO Q 

AW(v*)^]T(^cg fc A<(v*), A ( n V)^-D-^) fc T *— /o(v*), (80) 
fc=i k=a x 



MU-v*) = M* (v*) 



2 



T* ( V 
,* , J l I v 



(81) 
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and we have already specialized to Maxwell molecules, so that v* = p* /T*. In order 
to apply the consistency conditions (52) in the derivation of the hydrodynamic fields 
p*,u*i, and T-j*, it is convenient to define the moments 

$ = $(!) 4. <j>( n ) (ao) 

= / <*v* C ni V ;%r 3 A (I) (v*), (83) 

<ln 3 = / dv * ^o m »;" ! CA (II) (v'). (84) 
Therefore, the consistency conditions are 

$ooo = 0, (85) 

$ioo = 0, (86) 

$oio = 0, (87) 

$200 + $020 + $002 = 0. 



The evaluation of $nin 2 n 3 and $n 1 I n 2 n 3 is carried out in Appendix A. The first- 



order profiles are 



p* 1 (s\a)= P ^ ) (a)s, (89) 

u* !l (s\a)=ufl(a)s 2 , (90) 

T?(s\a)=Ti 3) (a)s 3 , (91) 

p{ 1 \a) = -Ul-^)T?\a), (92) 



where 

ni^ta) = - 

T &), , 4 a 4 7 F 2 (F\ + 2F 2 ) - Fi (1 - Fp) - 6F 2 

(a)= 3 a7Fo D(a) ' (94) 

with 

D(a) = 2 7 F 1 [F 2 -2F + Fi-27(Fo-2Fi)(Fi+2F 2 )]+a 2 FoFi(l-F ) 

-2a 2 7 [F (F 2 + 6FiF 2 + 8F 2 2 ) - 2F 2 (Fi + 2F 2 )] . (95) 

In the above equations the functions F r are understood to be evaluated at x = 7. 

As shown in Appendix A, the moment $ nin2 „ 3 is a polynomial function of s 
of degree n.i + n 2 + n 3 — 1. In particular, the non-zero elements of the first-order 
pressure tensor are 

P^Ma) = P^(a)s, (96) 

with 

^i(«)=3p^ 1) («)-Fi 1 ! 1 (a), (97) 
i£!iW = 0, (98) 
FE?i (a) = - + F 1+ 2F 2 ) T< 3) - [2 7 (fi + 2F 2 ) - 1] p« , (99) 

^i(«) = 1. (100) 

As for the heat flux vector, the results are 

qi 1 (a\s)=$l(a) + q %(a)s 2 , (101) 
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where 



1 

~2 

w 



2F X - F - 1 , F-L + 26F 
8^ 
1-Fq 
47 

Fx + a 2 



2 , 2 
a 



1 - F 7 Fx - 16F 2 



2 7 



,( 2 ) 



4 7 2 
47 



F 



(3) 



3F 2 



Fi — F 2 

7 



107 



(F X + 2F 2 ) 



(102) 



2J1 °'/ 2 ' 1 +F 1 +2F 2 -a 2 
2 7 

F+14F 2 



f9F -2Fx-7 



8 7 2 



4 



r, 



(3) 



-Fx - F + 2 7 (3Fi + 4F 2 ) 



-3a 2 



2 fl-Fh 1 
4 7 2 



Fx - 3F 2 



+2 7 (Fi + 2F 2 ) + a 2 



2 ^-^o , 1 
47 2 



Fi - 5F 2 



■J (4 + 4F + 5Fj + 2F 2 ) + ^ 7 2 (F X + 2F 2 ) + Sl. ( Fl - F ) . 
3 2 



(103) 



4%) 



1 - F 7 Fx + 2F 2 2 3Fx - F 2 - 2F 
3 



F, + Fo + a" 



4 7 2 7 2 
Fi — Fi \ (!) 2a 



T 



(3) 



Fi — F ( 2 ) 



2~ 



Pi 



-(2F 1 +F 2 ), 



(104) 



4 2 i(«) 



Fi - 1 1 a 2 
— 1 — Fi + F 2 H 

4 7 2 4 

-a (F0+F1 -2F 2 )4 2 { + 



1 - F 2F - 7 Fx + 14F 2 



2 7 2 



1 



(3) 



-(1-F )+ 7 (F 1 + 2F 2 ) 



-_-(F -3Fi + 2F 2 ) 



P W + £ [l - Fo - 2 7 (Fi + 2F 2 )] 
o 



(105) 



Equations (98) and (100) are consistent with the momentum balance equations 
(15) and (16), respectively. The energy balance equation (17) requires that 



ifl = a [ F o u x'i - 



(2) 



(106) 



Taking into account Eqs. (92)-(94), it is possible to check that Eqs. (105) and (106) 
are indeed equivalent. 

Let us now get the relationship between the scaled space variable s and the true 
(dimensionless) coordinate y*. From the definition (27) we have 

ay = T*( S \a,g*) 
ds p*(s\a,g*) 

= T*(s\a) + [Tx% S \a)-T*( S \a)pU S \a)]g* + 0(g* 2 ). (107) 
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Inserting Eqs. (60), (89), and (91) one gets 

y*(s\a,g*) = s - £ { 2 p«(a) - [if >(a) + 7 (a)p«(a)] s 2 } .9* + 0(g* 2 ). 

(108) 

4.2.2. Limit of small shear rates. As done in the case of g* = 0, it is illustrative to 
obtain the first-order coefficients (92)-(94), (97), (99), and (102)-(105) in the limit 
of small shear rates. Taking into account Eqs. (63) and (72) one gets 

rf ) («) = -^(l-fa 2 + ---), (109) 
1 29a 2 A 23136 2 



<i(a) = ---—- !- «' + ... , (110) 



According to Eqs. (40), (41), and (47), the NS equations yield v> x [ = -\, T[ 



2 5 V 725 
Tf) (a) = |(l + ^ + ...), (Ill) 

< ) i(«) = -f (l-^« 2 + -), (113) 
(0)/ , , 216a 2 / 23329 2 \ 

(2) , , 29a 2 / 1844 2 \ , 11C , 

(o), . 19a / 8172 , \ 

g- (^1 - "gg-a 3 + • ■ • J , (116) 

4%) = -a(l + 4a 2 + •••). (117) 

(2) _ 1 ^(3) 

anil i*±i ^, uiic cm-iciLiuiio _yit;iia (x 

^a (taking Pr = 1), and q^l = —a, so that they agree with the leading terms 
in Eqs. (110), (111), and (117). On the other hand, the NS approximation fails 
in accounting for the non-zero values of Pi , P^xii ^zz i> Q^i' QxV ano - % l- 1 R 
particular, q ®\ ^ even in the pure Poiseuille flow (a = 0) [40, 44, 55, 56]. 

5. Results and discussion. 

5.1. Finite shear rates. First order in g*. The nonlinear dependence on the 
reduced shear rate a of the zeroth-order quantities has been analyzed elsewhere 
[6, 19, 20], so that here we focus on the first-order corrections. Figure 2(a) shows 
the coefficients associated with the hydrodynamic profiles, i.e., (a), u^\(a), 

and Ti (a). The first two quantities arc negative, while the third one is positive, 
in agreement with what might be expected in view of Eqs. (109)— (111). On the 
other hand, the practical range of applicability of the truncated series (109)-(111) 
is restricted to small shear rates (a < 0.1). The addition of further terms in the 
(Chapman-Enskog) expansion in powers of a would not improve that range because 
of the asymptotic character of the series. Note that the range of applicability of the 
NS description (according to which p^ = 0, = — |, and = j^a) is even 
much more restrictive, especially in the case of p^ ■ 
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Figure 3. First-order coefficients q x l , q x [ (left panel), q y ll and 

Qy i (right panel) as functions of the reduced shear rate a. The 
dotted lines represent the terms shown in Eqs. (114)-(117). 



The normal-stress coefficients P^i( a ) an< ^ Pzl\( a ) are plotted in Fig. 2(b). Both 
coefficients vanish in the NS description. Again, the truncated series (112) and (113) 
are reliable only for a < 0.1. We observe that the xx element has a much larger 
magnitude than the zz element. The other relevant coefficients of the pressure 
tensor are not plotted because they are identically i = and Py y ]i = 1, as a 
consequence of the exact momentum balance equations (15) and (16). 
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Figure 4. Profiles of (a) temperature, (b) elements of the pressure 
tensor, and (c) flow velocity and components of the heat flux vector. 
The value of the reduced shear rate is a = 1. Two values of the 
external force are considered: g* — (dashed lines) and g* = 0.1 
(solid lines). 



The coefficients associated with the heat flux vector are plotted in Fig. 3. Ac- 
cording to the NS equations, q£\ = = = and = —a, what strongly 
contrasts with the nonlinear behavior observed in Fig. 3. This is especially dramatic 

(2) 

in the case of the strcamwise coefficient q], [, which deviates from zero even in the 

limit a — >• 0. It is interesting to note that, while q^\, q^\, and q^\ have definite 

signs (at least in the interval < a < 1), the coefficient q^\ changes from negative 
to positive around a ~ 0.42. 

In order to illustrate how the Couette-flow profiles are distorted by the action of 
the external force, we will take a = 1 with g* = and g* =0.1. While the latter 
value is possibly not small enough as to make the first-order calculations sufficient, 
our aim here is to highlight the trends to be expected when the fully nonlinear 
Couette flow coexists with the force-driven Poiseuille flow. The profiles are shown 
in Fig. 4. In the pure Couette flow (g* = 0) the temperature and pressure profiles 
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are symmetric, while the velocity and heat flux profiles are antisymmetric. The 
application of the external force breaks these symmetry features since the first- 
order terms have symmetry properties opposite to those of the zeroth-order terms, 
in agreement with the signs of S g in Table 1. As a consequence, the temperature 
gradient increases across the channel with respect to that of the Couette flow, as 
shown by Fig. 4(a). The elements of the pressure tensor are no longer uniform 
but exhibit negative gradients, especially in the case of the normal stress P* x [see 
Fig. 4(b)]. An exception is P* y , which is exactly uniform as a consequence of the 
momentum balance equation (15). To first order in g* the value of P* y is the same 
as in the Couette flow [see Eq. (98)] , but this situation changes when terms of order 
g* 2 are added [40, 55, 56]. We observe from Fig. 4(c) that the flow velocity (in the 
reference frame moving with the midplane y = 0) is decreased by the action of the 
external force. A similar behavior is presented by q*, what qualitatively correlates 
with the increase in the temperature gradient observed in Fig. 4(a). Regarding the 
component q*, it takes larger values in the Couette-Poiscuillc flow than in the pure 
Couette flow. In both cases (g* — and g* = 0.1) the shearing is so large (a = 1) 
that the two components of the heat flux have a similar magnitude, i.e., ~ \q*\. 
An interesting effect induced by the external field is the existence of a non-zero heat 
flux at y* = 0, even though the temperature gradient vanishes at that point. More 
specifically, q^(O) = q^\ > and q y {0) = q y °^ < 0. The sign of the former quantity 
changes, as noted before, at a ~ 0.42. 

5.2. Small shear rates. Second order in g*. Thus far, all the results are valid 
for arbitrary values of the reduced shear rate a but are restricted to first order in 
the reduced external field g*. One could continue the perturbation scheme devised 
in Sec. 4 to further orders in g* but the analysis becomes extremely cumbersome 
if one still wants to keep a arbitrary. Furthermore, the perturbation expansion in 
powers of g* is expected to be only asymptotic [56]. On the other hand, we can 
combine the results obtained here [see Eqs. (72)-(78) and (109)-(117)] with those 
of Ref. [56] to obtain the hydrodynamic and flux profiles to second order in g* , by 
assuming that both the reduced shear rate and the reduced force are of the same 
order, i.e., a ~ g*. The results are 



*/ I *\ i ^ * , u *2 2 , 

P (s\a,g ) = 1 - —ag s + -g s + 
u x {s\a,g ) = as- -g s H , 



(118) 



(119) 




(120) 



(121) 



Pyy(s\a 7 g*) = 1 - - 



a 



2 




(122) 




(123) 



Px V (s\a,9*) = 

?x(shs*) 



a + g*s + ■ ■ ■ , 
-9* + ■ ■ ■ , 



(124) 
(125) 




(126) 
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Figure 5. Profiles of the scaled temperature difference (T* — 
l)/g* 2 for \g*\ -c 1 and several values of the ratio a/g*: a/g* = 

(-), a/g* - 1 (- - -), a/g* = ^19/5 ( ), a/g* =3 (■■■), 

a/g* = 2^19/5 (- •■), and a/g* = 5 ( ). The left panel 

corresponds to the kinetic-theory predictions, Eq. (129), while the 
right panel corresponds to the NS predictions, Eq. (41). 



In these equations the ellipses denote terms that are at least of orders a 3 , g* 3 , a 2 g* 
or ag* 2 . As for the relationship between y* and s, from Eqs. (27) or (29) we get 

s = y* + -^a 2 y* 3 - ^ag*y* 2 (36 + y* 2 ) + ^<7* V* (22 + y* 2 ) +■■■. (127) 

Therefore, we can safely replace s by y* in Eqs. (118)-(126). 

Again, it is instructive to compare Eqs. (118)— (126) against the NS predictions 
worked out in Sec. 3 (with Pr = 1, in consistency with the BGK value of the 
Prandtl number). As can be seen from Eqs. (35), (40), (41), and (44)-(47), the NS 
expressions do not contain terms of order higher than a 2 , g* 2 , or ag* and therefore 
their comparison with the kinetic-theory results (118)-(126) is not biased. We see 
that only the NS results for u* and P* are supported by kinetic theory. This 
does not mean that Newton's law for the shear stress, Eq. (32) is satisfied, since 
t] = pjv and the hydrostatic pressure p is not actually uniform. As said before, 
the NS constitutive equations also fail to account for the existence of normal stress 
differences (typical non-Newtonian effects) as well as of a streamwise component 
of the heat flux (failure of Fourier's law). Perhaps the most interesting and subtle 
differences refer to the presence in Eqs. (120) and (126) of the extra terms s 2 
and — ^-ag*, respectively, which are absent in their NS counterparts, Eqs. (41) and 
(47). The extra term in q* implies that 9^(0) ^ 0, what represents a clear violation 
of Fourier's law (34). The term ^g* 2 s 2 present in the temperature field (120) has 
dramatic consequences on the curvature of the temperature profile at the midpoint 
y* = 0. From Eq. (120) we get 



dy 



*2 



2 2 38 „2 



d 2 T* 

= -~a l + —g*\ (128) 
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Therefore, while the NS temperature profile presents a maximum at y* = [see Eq. 
(43)], Eq. (128) shows that the profile has actually a local minimum if a 2 < ^§-g* 2 - 
To analyze this feature in more detail, let us rewrite Eq. (120) as 

T*-l l/»\ 2 , 2l 2 / a\ , 3i , 2 /19 1 „ 2 



g* 2 ~ 5 \g*J V + 15\^) V +V {25- m y J- ^ 

Figure 5(a) shows the scaled temperature difference (T* — l)/<?* 2 , as given by Eq. 
(129), for a/g* = 0, 1, ^19/5, 3, 2^19/5, and 5. In the case a/g* = (pure 
Poiseuille flow) the temperature profile has a minimum at y* = surrounded by 
two symmetric maxima at y* = ±2^/19/5. When a/g* ^ (Couette-Poiseuille 
flow) several possibilities arise. If < a/g* < VW5 

one still has a local minimum 
at y* = but the two maxima are no longer symmetric: the one with y* < moves 
to the center, while the one with y* > departs from it. This is represented by the 
case a/g* = 1 in Fig. 5(a). At a/g* = ^/l9/5, the left maximum and the central 
minimum merge to become an inflection point of zero slope. Next, in the range 
-J19/5 < a/g* < 2^/19/5 the temperature presents a local maximum at y* = 
followed by a minimum and an absolute maximum, both with y* > 0. This situation 
is illustrated by the case a/g* — 3 in Fig. 5(a). At a/g* = 2^/19/5 the minimum 
and maximum with y* > merge to create an inflection point of zero slope. Finally, 
if a/g* > y/19/5 [see case a/g* — 5 in Fig. 5(a)] only the central maximum remains 
and the profile becomes more and more symmetric as a/g* increases. In the limit 
a/g* — > 00 (or, equivalently, g* —> 0) one recovers the pure Coucttc flow. This rich 
phenomenology is absent in the case of the NS temperature profile, as shown in Fig. 
5(b). 

Given the physical interest of Eq. (120) or, equivalently, Eq. (129), it is convenient 
to rewrite it in real units. This yields 

T(y) = T(0)~A( d ^) 2 y> + ^( d -^) gy* 



2k(0) V dy J y=0 a 3 K (0) V dy 



m 2 n 2 (0) 
12 ?7 (0) K (0) y V ^^ r jfc|r(0) 



where Ct = |§ in the BGK model, while Ct — 1.0153 in the Boltzmann equation 
for Maxwell molecules [40, 44, 55]. Equation (130) still holds in the NS description, 
except that Ct = 0. 

6. Concluding remarks. In this paper we have studied the stationary Couette- 
Poiseuille flow of a dilute gas. As illustrated by Fig. 1(c), the gas is enclosed 
between two infinite parallel plates in relative motion (Couette flow) and at the 
same time the particles feel the action of a uniform longitudinal force (force-driven 
Poiseuille flow) along the same direction as the moving plates. Our main goal has 
been to assess the limitations of the NS description of the problem and highlight 
the importance of non-Newtonian properties. 

In order to get explicit results, the complicated Boltzmann collision operator has 
been replaced by the mathematically much simpler BGK model with a temperature- 
independent collision frequency (Maxwell molecules). The kinetic model has been 
solved to first order in the reduced force parameter g* for arbitrary values of the 
reduced shear rate a. Moreover, complementing these results with those obtained 
in previous works for the pure Poiseuille flow to second order in g*, we have been 
able to get the solution to second order in both a and g* . 
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Starting from the pure nonlinear Couette flow, we have studied the influence of a 
weak external force, as measured by the nonlinear shear-rate dependence of the nine 
coefficients p« (a), u^(a), if' (a), P«», Pj^a), <g(a), q^(a), q<$(a), and 

q y x {(a). These functions arc plotted in Figs. 2 and 3. A more intuitive picture on 
the distortion produced by the force on the Couette profiles for the hydrodynamic 
fields (p, u x , and T), the pressure tensor (P xx , P yy , P zz , and P xy ), and the heat 
flux (q x and q y ) is provided by Fig. 4. 

Complcmcntarily, we have obtained the quantities of interest [cf. Eqs. (118)— 
(126)] when the shear rate and the force are treated on the same footing, both to 
second order. This has allowed us to analyze [see Fig. 5(a)] how, by starting from 
the pure Poiseuille flow, the symmetric bimodal temperature profile is strongly 
distorted by the shearing until arriving at the symmetric parabola characteristic of 
the pure Couette flow. 

To put the paper in a proper context, it must be stressed that in our treatment we 
have focused on the nonlinear bulk properties of the Couette-Poiseuille flow, while 
ignoring linear effects like Knudsen boundary layers. The latter effects, however, 
may be relevant in applications like slow MEMS flows where the nonlinear effects 
considered here might be less important. 

Considering the great current interest in the force-driven Poiseuille flow as a 
playground to test hydrodynamic theories and theoretical approaches, we expect 
that the work presented here may contribute to motivate further studies, both 
theoretical and computational, on the Couette-Poiseuille flow. 
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Appendix A. Evaluation of <&nln 2 n 3 and &nin 2 n 3 - From Eqs. (80) and (83) one 
gets 

=EE ( k V)j^%^ k I ^^o ni -^" 2+fc+ ^r^* (v *), 

l=0k=9 e V 1 J [Ul l >- J 

(131) 

where 6q = 1 and 0i = for £ > 1. In Eq. (131) use has been made of the 
mathematical relations 

A( s )(-d s ) k B(s) = E (fj (~ds) k ~ £ { [d e s A(s)] B(s)} , (132) 

dXY 1 - ^$Ko ni ^- (133) 
Next, using the Maxwcllian integrals 

J dv*V* ni v* y n *v* z n3 MUv*) = K ni K n2 K n . d T* {n ^ +n *- 2), \ (134) 
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where K n = (n — 1)!! if n = even [with the convention (—1)!! = 1], being zero if 
n = odd, Eq. (131) becomes 

rtv(i) is \ ^ \ " fk + £\ ni\(-a) e k *( k+ni+n2+n 3 -2)/2 

k + n-i+n2 + n 3 - 2 T* 



=0 k=t 



Pi 



1 



(135) 

Before considering the integral ^nP,n 2 ,n 3 , it is convenient to make use of the 
relation 



(-d s ) k [A(s)B(s)] = (*) [(-d s ) m ^)] [(-5 s ) fc - l S( S )] (136) 

7Ti=0 ^ 7 

to rewrite the function A ( - II - l (v*) as 

= (1 - 7S 2 )A( n ' 0) (v*) + 2 7 sA( II ' 1) (v*) - 2 7 A (II < 2 >(v*), (137) 

where 



fc=0 



In the last step we have made use of Eq. (56) and the mathematical property 

£Cr)-pi:r> <-> 

e=o v 7 v 7 
Insertion of Eq. (137) into Eq. (84) gives 

$(") = fl - -vs 2 )^ 11 ' ) + 2<7s<J> (IU) - 2-/$ (II < 2) (1401 

where 

$ (II,m) = /" dv *y* «i w *«2 *n 3A (II,m)/ *n /j^-v 

Using again Eqs. (132)-(134), one gets 

( „, m) ^ fk + £ + m+l\fk + £\ ni \(-af 

e=o k=o v / \ / v i / 

x^ fc+ „ 2+ , +m (-9 s ) fc T* (fc+ " 1+ " 2+ " 3+m - 3) / 2 . (142) 

Once the integrals $nln 2 n 3 and $i I 1 I |i 2 „ 3 are expressed in terms of Tq (s), Pi(s), 
u* 1 (s), and T^s), we can apply the consistency conditions (85)-(88) to get the 
hydrodynamic profiles to first order in g* . To that end, we first guess the polynomial 
forms (89)-(91), so that only the coefficients p{ \ u^\, and T± remain to be 

determined. It is straightforward to check that $qq = <f>g Q = 0, and thus Eq. (85) 
is identically satisfied. The remaining relevant quantities in Eqs. (86)-(88) turn out 
to be 

*& = 2aT^^l + 2(4 2 i + a£>)F u (143) 
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(144) 

(145) 
(146) 

(147) 

(148) 

(149) 
(150) 

(151) 

(152) 
(153) 



In the above equations use has been made of Eqs. (63)-(65). 

Inserting Eqs. (144) and (149) into the consistency condition (87) one simply 
gets Eq. (92). Next, insertion of Eqs. (143) and (148) into Eq. (86) allows one to 
obtain Eq. (93). Finally, use of Eqs. (145)-(147) and (150)-(153) in Eq. (88) yields 
Eq. (94). Note that from Eqs. (92) and (146) one gets s" 1 ^ 1 



020 



Taking into account that Tq , p\, u* i, and Tj* arc polynomials in s of degrees 2, 

1, 2, and 3, respectively, Eqs. (135) and (142) show that <&n\n 2 n 3 and ^nl'rws are 
polynomials of degrees ri\ + n 2 + — 1 and n\ + ra 2 + 713 + m — 3, respectively. 
Consequently, the moments defined by Eq. (82) are polynomials of degree n\ +n 2 + 
n 3 - 1. 

Let us proceed now to the evaluation of the integrals <f>n}n 2 n 3 and <&nln 2 n 3 related 
to the pressure tensor and the heat flux. The integrals related to the diagonal 
elements of the pressure tensor are already given by Eqs. (145)-(147) and (150)- 
(153). For instance, Py V .\ = P* + $020, with similar relations for P* xl and P* zl . 
The results are displayed in Eqs. (96)-(99). In the case of the shear stress, one has 
P* j = $ 110 . From Eqs. (135) and (142) we obtain 



v no 



T 



(11,0) 

110 



(3) „ Fp — 3Fj + 2F 2 

7 

(11,0) 



2 Mi 2 iFi -p\ ± 'a(2F 1 -F ), 
2 



„(i). 



0, $ 



110 



F, $ 



(11,2) 
110 



;(F 1 -F 2 )S. 



(154) 
(155) 



This gives 



F* Ms) 



F - 3F + 2F 
1 

7 

§ 7 (Fi+2F 2 ) ; 



T 



(3) 



a{2F 1 - F )p\ 



(i) 



(156) 
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Making use of Eqs. (92) and (93), it is easy to check that Eq. (156) reduces to Eq. 
(100). 

In the case of the heat flux, one has 

<i = 5 ( $ 2i" + $ "3o + $012) - ^, <i- (158) 
After tedious algebra one obtains Eqs. (101)-(105). 
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